Hartree-Fock-LAPW method: using the full potential treatment for exchange 
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We formulate a Hartree-Fock-LAPW method for electronic band structure calculations. The method 
is based on the Hartree-Fock-Roothaan approach for solids with extended electron states and closed 
core shells where the basis functions of itinerant electrons are linear augmented plane waves. All 
interactions within the restricted Hartree-Fock approach are analyzed and in principle can be taken 
into account. In particular, we have obtained the matrix elements for the exchange interactions of 
extended states and the crystal electric field effects. In order to calculate the matrix elements of 
exchange for extended states we first introduce an auxiliary potential and then integrate it with an 
effective charge density corresponding to the electron exchange transition under consideration. The 
problem of finding the auxiliary potential is solved by using the strategy of the full potential LAPW 
approach, which is based on the general solution of periodic Poisson's equation. Here we use an 
original technique for the general solution of periodic Poisson's equation and multipole expansions 
of electron densities. We apply the technique to obtain periodic potentials of the face centered cubic 
lattice and discuss its accuracy and convergence in comparison with other methods. 
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I. INTRODUCTION 

In contrast to successful applications of the Hartree- 
Fock (HF) method to molecules, which have become a 
routine procedure, Hartree-Fock calculations of electron 
band structure of solids are relatively rare. The electronic 
band structure calculations mainly rely on the density 
functional approachau with the local density approxima=j 
tion (LDA) or other non-local exchange approximationsu 
which try to mimic the Fock exchange for conduction 
electrons. As known, LD approximation is a product 
of the Hartree-Fock theory of the free electron gas and 
as such it inherits its shortcomings .□ The LD approx- 
imation being universal and simple from one side, is 
structureless and leads to an unphysical self-interaction,El 
from the other. In addition, LDA is an uncontrolled 
approximation, and there is no straightforward path to 
further improvements in its accuracyfl Therefore, it is 
highly desirable to combine the Hartree-Fock approach 
with methods based on plane wave basis sets which tra- 
ditionally play an important role in the field of electron 
band structure calculations. Although a variety of com- 
putational procedures on xiexiodic HF method has been 
reported in the literature,ElEj all of them use Gaussian- 
type orbitals (GTO). In this paper we present a merger 
of the HartpepEock approach with the linear augmented 
plane waveoii3 (LAPW) method, which constitutes a 
new Hartree-Fock-LAPW procedure for electron band 
structure calculations. The consideration of the LAPW 
basis is important since the method has been extensively 
tested and it is traditionally associated with electronic 
band structure calculations. 

At the center of the Hartree-Fock-Roothaan approach 
there are calculations of the direct Coulomb and the ex- 



change matrix elementsEzI and in the present work we con- 
sider thoroughly all cases, including the delocalized and 
localized states. The advantage is that all expressions 
for the matrix elements can be found in one source (i.e. 
this paper), but on the other hand, this has led to nu- 
merous formulas and their extensive derivations. While 
GTOs are popular in calculations of molecular systems 
and have the advantage of giving easily integrable poly- 
ccntric functions, we lose this benefit working with the 
LAPW basis set. There, we have to rely on a different 
strategy of calculations. Our choice rmr^Js the full poten- 
tial LAPW (FP-LAPW) treatment j&B which is based 
on the procedure of finding the general solution for a pe- 
riodic Poisson's equation. For Poisson's equation we use 
a new o rigin al technique (Sec. II). It employs the Ewald 
methodE2lO for the potential of the monopole terms and 
an exact expansion in Fourier series in the interstitial re- 
gion, with a multipole series expansion of the potential 
inside the spheres. 

A principal new element of our work is to use the strat- 
egy of the FP-LAPW method for the calculation of ex- 
change matrix elements between itinerant electrons. In- 
deed, an exchange matrix element {ab\V ul \ba) (a and b 
here represent electronic states with the same spin com- 
ponent) can be estimated as a Coulomb self-interaction 
energy of the charge density distribution ip%(R)ip a (R), 
which appears as a result of the electron transition a —> b. 
The procedure implies (i) an introduction of an auxiliary 
potential of the charge density ipl(R)i/j a (R) and (ii) the 
integration of the potential with the charge density. Such 
treatment is rigorous and ensures that all contributions 
including the long range ones (like the exchange Ewald 
sums) are carefully considered. It is worth to notice that 
a simple truncation of the exchange series of other meth- 
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ods excludes such long range effects. 

Here we also consider crystal electric field (CEF) effects 
for the closed shell core electrons. Usually in the full 
potential LAPW treatment the core electrons are taken 
in the spherical approximation and the CEF effects are 
ignored.Ej While the crystal field splittings are small it 
may be still necessary to include these effects for precise 
calculations of solids where one compares the stability of 
different crystalline phases where the energy difference 
is of the order of fOO K. The present consideration's a 
generalization of previous works on CEF effectsr3c3 

For clearness and simplicity in the following we con- 
sider only one atom per unit cell. However, the method 
can be easily generalized for the case of a few atoms. In 
this paper we consider the core states being completely 
confined inside the spheres. Various improvements of the 
method in this respect are also possibleJIS but they are 
not of fundamental importance and will not be concerned 
here. 

The paper is organized as follows. We start with the 
general solution of periodic Poisson's equation, Sec. II. In 
Sec. Ill we illustrate how this technique works for the case 
of a face centered cubic (fee) crystal. We compare the 
convergence and accuracy of different approaches. Then 
we discuss the Hartree-Fock-Roothaan method for solids 
with extended states, Sec. IV, and calculate the direct 
Coulomb matrix elements for LAPW basis functions. In 
Sec. V we present calculations of the matrix elements of 
exchange, while Sec. VI is reserved for our conclusions. 



II. GENERAL SOLUTION OF PERIODIC 
POISSON'S EQUATION 

The general solution of Poisson's equation V Coul for 
a periodic charge density was demanded by necessity to 
improve-the-l'muffin-tin'' (MT) approximation of LAPW 
method.EfrHl LAPW electron band structure calcula- 
tion procedure with the general potential V tot (R) = 
V Coul (R) + V exc (R) is referrecL-to as the full poten- 
tial LAPW (FP-LAPW) methodic since in principle no 
shape approximation (except LDA) for the potential is 
assumed. First full potential calculations have been re- 
ported in Refs. ^,^6[ and since then they are often used 
to investigate the electron band structure of solids. These 
FPLAPW calculations are based on the general solution 
of Poisson's equation given by Weinert in Ref. [1^. The 
basic idea which has been proposed already in Ref. |l9| in- 
cludes two steps: I) to obtain the potential in the inter- 
stices and then 2) to solve the boundary value problem in- 
side a sphere. The point I ) ia.Rjef. p"8| represents an alter- 
native to the Ewald methodclJ'EJ where the true charge in- 
side the spheres is replaced with a pseudo-charge-density 
of the same multipole moments. In the present paper we 
propose a new technique which uses the Ewald expansion 
for the monopoles (I = 0) and the exact Fourier expan- 
sions for the higher multipoles (/ > 0) in the interstitial 




FIG. 1. The region inside the spheres S (white) and the 
interstitial region / (gray). 



region. The procedure is straightforward and besides a 
truncation of the multipole series, it contains only one 
cut-off parameter for the Ewald expansion. The s trate gy 
consists of the same two steps as outlined above.E3E3 



A. Dual and equivalent representations of charge 
density 

We start with the description of charge density in a 
crystal. For a periodic crystalline structure we introduce 
n to label the sites and the direct lattice vector X(n) 
which specifics the centers of atoms (spheres) . The posi- 
tion vector R inside a sphere S(ft) is given by 



R(n) = X(n) + r(n). 



(2.1) 



In spherical coordinates r(n) — (r(n), f (n)), where f = 
fl = (0,(f)). For a Bravais lattice of X(n) we construct 
the reciprocal lattice of vectors K. It is convenient to 
partition the space into the region S inside the non- 
overlapping spheres of radius R and the interstitial re- 
gion /, Fig. 1. We expand the charge density in Fourier 
series in region / and, in multipolar terms in region S 
(dual representationEaEj) : 



Rei * — ' 



iKR 



ep(R) 



Res 



(2.2a) 



K 



E ?A(r) S A (f) - Z ^=^SoJ • (2.2b) 



Angular functions S\ are symmetry adapted functionsE£l 
(SAFs) which are linear combinations of spherical har- 
monics (see Appendix A for details and definitions). The 
total density is an invariant of the space group under con- 
sideration. Therefore, in the expansion (2.2b) we use only 
SAFs Sa of the full (or unit) A\ g symmetry. Finally, 5(r) 
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in fl2.2b|) is the one dimensional (radial) delta- function, Z 
is the atomic number and e < is the charge of electron. 
In the following we will use atomic units and the charge 
density i s und erstood in units e = — 1. The densitjj-efc 
pansion (2.2a,b) called also the dual repj-esentationE3'Ea 
is characteristic of the LAPW method. E£li!£l In fact, the 
dual representation is natural because in the interstitial 
region electron density is smooth whereas large charge 
oscillations near nuclei require a multipole analysis. 

Using the linearity of Poisson's equation we can extend 
smoothly the interstitial Fourier charge representation in- 
side the spheres and subtract it out. We obtain 

p(R)=p I (R) + '£p' s(H) (R), (2.3a) 
Ps(n)(R) - (ps(n)(R) ~ Pi(R)) ©(£ e S(n)), (2.3b) 



where pi{R) is defined by ( 2.2a ) over the whole space, 
P's(n) ^ S t Qe renormalized charge density inside a sphere 
S(rt) and 9 is the unit step function. Using Eq. ( |2.l| ) and 
the plane wave expansion in SAFs centered at the site n, 

el KR = e iXX(X) An J2i l Ji(Kr(n))S A (K)S A (r(n)), (2.4) 

A 

we rewrite the multipolar density inside the sphere as 



Po(r) = Po(r) - VAttpi(K = 0)-Z 



Sir) 
Aur 2 



(2.5a) 



PA(r) = P\(r) - 4iri l J2 jl(Kr)S A (K) Pl (K). (2.5b) 



We s hall refer to (2.3a,b) with the multipole densities 
( 2.5a ,b) as to the equivalent representation of the charge 
density. We will use it every time when we need to obtain 
the pot entia l in the interstitial region because the Fourier 
series (2.2a) is now extended to the whole crystal, that 
has some important advantages (see next subsection). 



In principle, in Eq. ( 2.2a ) there is a term with K = 
that corresponds to the homogeneous electron density 
distribution. However, it will not contribute to the re- 
sulting potential because it cancels with the positive ho- 
mogeneous density distribution arising from nuclei (see 
also next subsection). 



B. The solution in the interstitial region 

The potential in the interstices has three components. 
The first contributio n aris es from the density pi (R) given 
by the Fourier series (2.2a) throughout the whole crystal. 
The second and third terms are due to the monopole (I = 
0) and the high multipole (I > 1) moments of electron 
density inside the spheres. (We recall that here we are 
dealin g with the equivalent charge density given by Eqs. 
~5gb).) 
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From the Fourier series of pi(R) we find 
Vi(R) = ^V(l?)e iJ?i? , 

K=£0 



where 



V\K) = ^pj(K). 



(2.6) 



(2.7) 



The easiness of the last expression owes to the fact that 
we extend the Fourier expansion of pi(R) to the whole 
crystal. A complication is that the density inside a sphere 
is renormalized according to Eqs. (2.5a,b). Next we con- 
sider a sphere S(n) and for a moment we will not include 
to the potential the contribution due to the charges out- 
side the sphere. Thus we introduce a single sphere potep=. 
tial Us{r) = Us(n)(r(n)). Using the one-site expansionE2l 
in SAFs 
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\r(n) — r'{n)\ 



= v A\(n = fj', r, r')SA(r')SA(r), (2. 



where 



VAA' (n = n' , r, r') — / dfl / <if2 



v Sa(Q)Sa>(SV) 
\r(n) — r 1 (n)\ 



= Saa> 



4-7T r' K 
21 + 1 rl +1 : 



(2.9) 



and r< (r>) is the smaller (larger) of r and r' , we obtain 
the solution inside the sphere (0 < r < R), 



Usif) 



Mr) 



Yo(r) 



A#0 



Here we have introduced the functions 
47rr 3 ^ 

Qo(r) = qo(r) - —Pi( K = 0). ( 2 - lla ) 

Q'o(r) = <l'o(r) - MR 2 - r 2 ) Pi(K = 0), (2.11b) 
with 

pr 

(2.12a) 

q' (r)=Vbr I p (r')r'dr', (2.12b) 

J r 

for the spherically symmetric term A = I = Q, and 

(2.13a) 
(2.13b) 



qo (r)=ViTr / p Q (r')r"dr' - Z, 
Jo 

r R 



dr , 
1 1-; 



Qa(0 = / p A (r )r 
Jo 

Q A (r)= l R p' A {r')r' L - l dr', 
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for the other multipoles, A / 0. Ta king into account 
t he ex plicit expression for p'(r), Eq. ( 2.5b ), we rewrite 
fl2.13a| ,p) as 



Q A (r) = q A (r) - 47r i y+ 2 'Sa{K) Pi {K) 



jj+iCgfO 
K 



Q' A (r)^q' A (r)+4m l J2 S A (K) P i(K) 
1 (h-xiKR) ji^(Kr 



where 



K 



9a (r) 



R l - 



-1-1 



p A (r)r dr , 



Q A (r) = I PA(r')r fl - l dr' 



(2.14a) 

(2.14b) 

(2.14c) 
(2.14d) 



In order to obtain Eq. ( [2.14a ,b) we have performed inte- 
gration of spherical Besscl functions using the properties 



10.1.23, 10.1.24 of Ref. 31 



For the potential Us(r) outside the sphere S(n) (r > 
R) we have 



Qo(R) 



^ 21 + 1 H+i bA[r) - 



A#0 



(2.15) 



Now we are ready to compute the potential of the peri- 
odic arrangements of the spheres, 



V s {R)=J2Us(R-X(n)). 



(2.16) 



Since we seek for the solution in the interstitial region it 
is desirable to expand it in a Fourier series, in the spirit 
of the dual representation, 



V S (R) = ^'v s {K)e liiii . 



(2.17a) 



We distinguish two contributions, Vq and Vfj, from the 
spherically symmetric charge distribution and from the 
other higher charge multipoles, respectively: 



V s (K) = V s (K) + Vi I (K). 



(2.17b) 



The monopole term (I — 0) Vq can not be computed 
directly and requires a seeeial treatment. Here we em- 
ploy the Ewald methodElO and obtain for the Fourier 
coefficients 



v K A 



MR), 



(2.18) 



where G is the Ewald cut-off parameter and v is the unit 
cell volume, v = V/N. The improved Ewald technique re- 
quires an additional summation in real space, see Ref. \2~U. 



Since we are concerned with the Fourier expansion in the 
interstitial region we can always take G large enough and 
neglect the summation in real space. (At large G the 
dimensionless parameter c = GR 3> 1 and the comple- 
mentary error function decreases exponentially.) Then 
the Ewald expansion with the Fourier coefficients ( 2.1S| ) 
converges absolutely in the inte rstiti al region. 

In the Ewald expa nsion, Eq. ( 2. IS ), and in the Fourier 
decomposition (2.6) we have omitted the terms with 
K = which correspond to a homogeneous charge distri- 
bution of the electrons and the nuclei. Both these terms, 
considered separately, diverge. Their sum, however, as a 
consequence of electroneutrality gives the zero value of 
the charge density and therefore can be discarded. 

The high multipole terms (/ > 1) converge fast. Their 
Fourier coefficients are given by 



V&8) = \l e 

v ./cell 



-iKR 



V SM (R)d 3 R. 



(2.19) 



Here the integration is over any primitive unit cell of the 



crystal, and Vs M (R) is the potential (2.16) where only 
the high multipoles ( I > 1 ) are retained in Us- Since the 
result of integration ( J2.19 ) does not depend on the choice 
of unit cell, it followsEi 
be rewritten as 



that the coefficients Vf t (K) can 



V m(K) = - 



~iKR 



V J crystal 



U SM (R)d 3 R 



(2.20) 



Here the integration is taken over the whole crystal and 
Us M {R) is given by the last (high multipole) term on the 
right hand sides of ( |2.10[ ) and ( |2.15| ) for < \R\ < R and 
\R\ > R, respectively. Using Eq. 10.1.24 of Ref. |l] we 
find 



viAK) 



where 



(47T) 



-v- 

A^O 



-S A (K) 



x [A A (K) + Q A (R) 



KR- 1 J 



A A (K) = J (90- + r l Q' A (r)^j 3l (K r ydr. (2.21b) 



(Note that at r — > 0, Q/ 



J+2 



ji(Kr) 



and the integral ( p. 21b ) does not diverge.) 
Co lle cting all contributions together, Eqs. ([2.7|), 
( |2.18 ), (2.21a), we obtain for the Coulomb potential in 
the interstitial region, 



V(Rel) = J2 v (*)< 



iKR 



where 



V(K) = V Z (K) + V S (K) + V&(K). 



(2.22a) 



(2.22b) 
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C. The solution inside the spheres 



The potential in the form (2.22a,b) can not be ex- 
tended directly to the region S inside the spheres. The 
main reason is that it requires an infinite number of plane 
waves to describe the Coulomb nucleus singularity. In 
order to obtain a well-behaved solution we expand the 
total potential yCSl-inside a sphere S(n) (r(n) < R) in 
a multipole serieslijll3 



V(f(n)) = V (r) + 5^V A (r) S A (f), 

A#0 



(2.23) 



The potential V{r(n)) is defined uniquely by boundary 
conditions.EJ These are of the Dirichlet kind here, since 



by using the expansion ( 2.22a ) for V(R S I) on the 
sphere surface, we arrive at 



In the following we drop the argument R in the boundary 
values V = V {R) and V A = V A (R). Notice that at r = 
the potential due to the charges outside the sphere S(n) 
is given only by V out since the other contributions (I ^ 0) 
reduce to zero. 

Recalling t hat the Fourier components of the intersti- 
tial potential (2.22a,b) consists of three parts, we distin- 
guish the same three contributions on the sphere surface 
for V and V A in Eqs. ( |2.24a| ,b): 



Va = v! + v, 



So 



V? 



(2.29) 



Here V£, V£° and V£ M are obtained by using Eqs 



f° and Vl M 

( |2.24a| ,b) where instead of V(K) we substitute V T (K), 
V n s {K) and V^(K), respectively. Taking into account 
( |2.28a ,b) we can separate the contributions to V A ut into 
three parts, 



V (R) = J2'i°(KR)V(K), (2.24a) 
V A {R) = 4Tri l Y,'ji(KR)S A {k)V(K). (2.24b) 

Now inside-, the sphere S(n) we have a Dirichlet 
problem ESEj We solve it by-, considering the spherical 
Green's function expansion,E2l and rewrite the solution 
(Eqs. (3.126) and (3.114) of Ref.|o|) in the following form: 



V(f(n)) = U s (r(n)) + Vg ut (r(H)), 



(2.25) 



The potential Us{r(n)) is induced by the charge distri- 
bution inside the sphere S(n). Clearly, it is independent 
of site n and we find 



Qa{r) 



A/0 



47r / q A (r) 



21 + 1 



+ r l q' A (r))S A (f). (2.26) 
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Not ice th at the potential U s diffe rs horn Us, Eq. (p.l5[) . 
In (|2.26D we use q A , Eq. (fe.l4cj ), and q' A , Eq. fl2.14d| ) 
(not Q A and Q' A \)- q A and g A corres pond to the initial 
charge density inside the sph ere, Eq. ( 2.2b| ), rather than 
the equivalent one, Eqs. ( 2.5a , b). Vg ut (r) is the potential 
due to the charge situated outside the sphere S(n), 



TrOUt 

V S 



(r)=^(i?) + ^V / 

A/0 



out 
A 



r V 



(R)[-) S A (f). (2.27) 



Here the constants V out 



are given by 

V'OUt 




V'OUt 
A 



V - 

v A 



V out (R) and V£ l 

qo(R) 

R ' 

4tt q A (R) 



21 + 1 



V£ ut (R) 



(2.28a) 
(2.28b) 



V'OUt TrOUt 
A — V AJ 



V, 



A, So 



v; 



out 

A,Sm' 



(2.30) 



Here V A U / is the contribution from the interstitial re- 
gion, while Vj^g and V^"§ accounts for the contribu- 
tions from the monopoles and the higher multipoles from 
the sites n' ^ n, respectively. Explicitly, we obtain 



(2.31a) 
(2.31b) 

(2.31c) 



\rout 
V A,I " 


-vl, 






\rOUt 

V A,S 


= (V S °- 


qo(R) 

R ' 


A = 0, 
A^O, 


V A,S M 


r v s ", 




A = 0, 


\ yl M 


4-7T 

21 + 1 


ri+i > a r u - 



The representation of V(R) given by Eqs. ( 2.25| ) and 
(|2.30|) is convenient because one clearly distinguishes the 
intra-site part Vs from the inter- site ones, Vjfg , V^g , 
and Vgf . 



D. Comparison with the two-center expansion of the 
Coulomb potential 

It is instructive to compare the results of the previous 
subsectio n wit h the two-center expansion of the Coulomb 
potential,BEl 

1 



\R(n) - R'{n')\ 

= $>aaW; r, r') S A {f{n)) S K > (r'(n')), (2.32) 



AA' 



where 

r- >\ fjnf~\ (jcif~'\ S A (n)S A ,{h') 

v AA i(n,n ; r,r ) = dll(n) I d\l{n ) 



\R{n)- R'(n')\ 

(2.33) 
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The potential o n the sphere surface at a site n is found 
by integrating ( [2.32 ) over r'(n') and summing over the 
spheres n'(^ n): 

V s ° ut (Q) = V° ut + EVa,1 + Vg£ M ) S A (Q), (2.34) 

A^O 



where (compare with Eqs. (2.31b,c)) 



V'OUt 
A,Si 



v A0 (n,n',R), 
V$ M = E'E' w AA'(n,n',E), 



(2.35a) 
(2.35b) 



and 

/•« 

v AA >(n,n' ,R) = / v AA/ (n,ri'; R,r')p' A ,(r'y 2 dr' . (2.36) 
Jo 

Sincei 

v AA '(n,n'; r,r') 



(r)'(r') 



|X(n) -X(n')| i+i ' +1 



(2.37) 



the sums on n', Eqs. (2.35a,b), converge fast (except for 
the case I = I' = which is irrelevant here). Practically, 
it is enough to consider summations over a few neighbor- 
ing shells. 

In many cases (f or the cubic symmetry, for exam- 
ple) the whole sum (2.35b) is small in comparison with 
( p.35a| ). Eq. ( 2.35a| ) then corresponds to the potential 
calculated with spherical symmetric charge distributions 
for the neighbors (i.e. A' = I' = 0) and a ho moge - 
neous charge distribution in t he interstices, Eqs. (2.5a), 
( 2.11a ). It follows from fl2.37 ) that v A o(n,n'; r,r ) does 
not depend on r' which is a dummy argument and 



v A0 (H,H',R) = B A (H')Q (R) 



4ir 



(2.38) 



Here Qn (R) is the effective charge inside the sphere given 
by Eq. (|2.11a ), and B A (n') — v A p( n, ft'; R,r'). Com- 
bining ( |2.35aD and (fp|) with fl2.27| ), we introduce the 
potential 



V CEF (r) 
where 



Qo(R)- 



A^O 



B A = Y^B A (ft'). 



(2.39a) 



(2.39b) 



The potential V represents the crystal electric field 
potential (CEF) for localized electrons at site n if there 
are no multipolar couplings with conduction electrons 
and the electron density in the interstitial region is ho- 
mogeneous, V (K) = 0. Examples of calculations of B A 



for a face centered cubic lattice (fee) will be given in 
the next section. In the nearest neighbor approxima- 
tion one has B A = N nn B A (ri'), where ft' refers to any of 
the nearest neighbors and N nn is the number of nearest 
neighbors. Such CEF potential has been considered in 
Refs. ^2 24. A more sophisticated CEF for localized and 
itinerant electrons will be discussed in Sec. IV. 



III. APPLICATION TO THE FACE CENTERED 
CUBIC LATTICE 

In this section we apply the technique of Sec. II to test 
the method and to calculate some important structural 
constants. As a Bravais lattice we take the face centered 
cubic (fee) lattice. 



A. Periodic Coulomb potential of the monopole 
(/ = 0) charge density 

We first consider the potential of the unit (q = 1) 
point charges situated on a face centered cubic lattice. 
For simplicity we take the cubic lattice constant a = 1. 
(The volume of the unit cell is v = 1/4.) We consider 
touching spheres on the fee lattice with the close contact 
radius R = \/2/4. As we know from Sec. II, the crystal 
potential is defined only if the total charge of the unit 
cell is zero. Therefore, complementary to the positive 
point charges we must introduce a compensating negative 
charge density distributed homogeneously in the intersti- 
tial region. The total charge confined in the interstitial 
region of the unit cell must be qi — — 1 with the density 
P/ = -1M = -4/(1-tt/3V2). 

Using the method of Sec. II we first expand the poten- 
tial of the interstitial region in Fourier series and then 
extend it smoothly to the whole crystal, Eqs. (2.17a,b). 
There are no high multipole moments inside a sphere 



and from Eq. (|2.21a| ) we find V§(K) = 0. Due to 



electroneutrality the terms with K = are omitted. 
Since in the interstitial region there is only one com- 
ponent pi(K = 0) which is different from zero, this 
results in V 1 (K 0) = 0. Therefo re, the potential 
in the interstices is given by ( 2.22a| ) where the only 
nonzero contribution to the Fourier component V(K) is 
the Ewald contribution V S {K), Eq. fl2.18| ). Notice that 
for the Ewald expansion we use an effective point charge 
Q = Q eff = q- Pl 4:TrR 3 /3 = 1 + tt/(3\/2 - tt). The 
renormalization of charges is the price that we have to 
pay for the extension of the Fourier expansion inside the 
spheres, Sec. II. 

The potential inside a sphere S(n) then is given by 



V(r) = - 
r 



j rout 



-K)4(^)V 4 (r) + U 06 (^)V 6 (r) 
+ .... (3.1) 



G 



TABLE I. Comparison of the Ewald expansion and the 
method of Ref. 18. G is the Ewald cut-off parameter, n is 
the parameter of Ref. 18. "Exact" corresponds to the Ewald 
expansion with K R < 1779 and G = 780, with the accuracy 
of » 10" 5 . 







Ewald 






Ref. |l| 




K-max R 


G 


Vo 


(Vi) 


n 


Vo 


(Vi) 


11.54 


7 


-0.4545 


-0.5607 


5 


-0.4996 


-0.5607 


14.74 


8 


-0.5131 


-0.6200 


8 


-0.5597 


-0.6200 


28.45 


12.5 


-0.6288 


-0.7359 


17 


-0.6287 


-0.7359 


46.33 


20 


-0.6778 


-0.7849 


39 


-0.6714 


-0.7849 


68.54 


30 


-0.6953 


-0.8024 


62 


-0.6849 


-0.8024 


90.73 


40 


-0.7014 


-0.8085 


84 


-0.6911 


-0.8085 


exact 




-0.70922 


-0.81630 




-0.70922 


-0.81630 



Here we introduce the cubic harmonics K& (O) and 
Kq (f2) which acaJSAFs S;=4,Aig(0) and Si = Q : Ai g (Q), 
correspondingly jESlj and 



V ° ut =Vo-^- 



(3.2) 



It is interesting to remark that although we have started 
with the model of point charges on the fee lattice, the 
total crystal potential (|]l]) acquires contributions of the 
high multipole orders, / = 4,6,8,.... The constants Vo, 
V04, Vo6 of the potential are found from the boundary 
conditions on the sphere surface, Eqs. (2.24a,b). The 
potential energy of the charges (fee lattice) is given by 

Ef cc /N = V out + -(Vr) = -3.9885 a.u., (3.3) 

where (Vj) is the average potential of the interstitial re- 
gion /, 

(Vi) = — f V(R)dR = — V 'v(K)(D(K). (3.4) 
vi J i vi ' 

The overlap integral O(K) reads 

4nR 2 ji(KR) 



0{K) 



K 



(3.5) 



of the 



Thus Ef cc /N, Eq. (3.3), has the meaning 
Madelung constant for the fee lattice. 

The values Vo, V04, Vo6 and (Vj) are proportional to the 
effective charge Q e f / ^ 1 ■ It is convenient to consider the 
normalized to the unit charge quantities Vo = Vo/Q e ff, 
Vq 4 = V M /Q e ff, V 06 = Vw/Qefj and (V £ ) = {Vi)/Q eff . 
We study the convergence of V*o, Vq4, Voe and (Vi) by 
employing the Ewald method and the technique of Wein- 
ert, Ref. |l8|. The results are quoted in Tables I and II 
where the same number of reciprocal vectors \K\ < K max 
is used for both approaches. (Though we did not op- 
timize the convergence of Fourier series by choosing the 
best G and n for the Ewald and Weinert approaches.) 



TABLE II. Comparison of the Ewald expansion and the 
method of Ref. 18. G is the Ewald cut-off parameter, n is the 
parameter of Ref. 18. 







Ewald 






Ref. y 




i^vnax R 


G 


v 04 


Voe 


n 


v 04 


Voe 


11.54 


7 


-0.18135 


-0.14304 


5 


-0.18233 


-0.14600 


14.74 


8 


-0.18186 


-0.14459 


8 


-0.18192 


-0.14459 


28.45 


12.5 


-0.18193 


-0.14467 


17 


-0.18192 


-0.14466 


46.33 


20 


-0.18192 


-0.14466 


39 


-0.18192 


-0.14466 


exact 




-0.18192 


-0.14466 




-0.18192 


-0.14466 



TABLE III. Calculation of V04 by summation over neigh- 
boring shells, N s h is the number of neighbors in the shell, 
n' stands for the coordinates of a shell representative, «04 is 
the two-center integral (2.33) for So and K4. "Exact" means 
summation over 128428 neighbors (d < 12) on the fee lattice. 



shell 


N sh 


m 


«04 


V04 


1 


12 


(1,0,0) 

(- - l) 

\ 2 ' 2 ' / 

(1,1,0) 
(!,io) 


-0.070694 


-0.23930 


2 


6 


0.049988 


-0.15470 


3 


24 


-0.004535 


-0.18540 


4 


12 


-0.002209 


-0.19288 


5 


24 


0.002782 


-0.17405 


10 








-0.18073 


20 








-0.18283 


exact 








-0.18191 



We observe that V06 and V04 converge faster than Vq and 
(Vi). 

The Ewald method estimates better the monopole 
terms Vq and (Vi), Table I. In this case it is simple and 
stable, while Weinert' procedure requires evaluation of 
spherical Bessel functions of higher order for each K. 

Both methods give almost the same results for V04 and 
V06, Table II. V04 and V06 can be calculated indepen- 
dently by means of the two-center expansion of Coulomb 
interaction, which we have considered in Sec. II D. In 
that case we first calculate VQ4{n,n') and voe(n,n'), Eq. 
( |2.33D , where r = r' = R, with 5 = l/V&r on a hrst site 
n, and the cubic harmonic (SAF) K4 or Kq on a second 
site ft'. Since the integrals vo4,(n,n') and voe(n, ft') de- 
crease very fast with the distance d between the two sites 
(l/d 5 for K4 and 1/d 7 for K§), we obtain V04 and V"o6 by 
summing vo4(n,n') and Voe(ft,n') over all neighbors ft' 
on the fee lattice, i.e. 



1 v-V 
Vn = —j=>^ v t{ft,ft'), 



(3.6) 



where t = 4 or 6. The results of such calculations are 
shown in Tables III and IV. The convergent values V04 
and V06 perfectly match those obtained by the techniques 
of Ewald and Weinert, Table II. However, the Ewald and 
Weinert' approaches require less computer time and are 
more efficient than the two-center expansion. Finally, 
we would like to remark that the calculated in Table II 
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TABLE IV. Calculation of V06 by summation over neigh- 
boring shells, N a h is the number of neighbors in the shell, 
ft' stands for the coordinates of a shell representative, U06 is 
the two-center integral (2.33) for So and K§. "Exact" means 
summation over 111956 neighbors (d < 9) on the fee lattice. 



shell 


N sh 


ft' 


"06 




1 


12 


(s,§,0) 
(1,0,0) 
(- 1 1) 

\ 2 ' 2 ' / 

(1,1,0) 
(1,1,0) 


-0.044247 


-0.14978 


2 


6 


0.002407 


-0.14571 


3 


24 


0.000299 


-0.14368 


4 


12 


-0.000346 


-0.14485 


5 


24 


0.000005 


-0.14482 


10 








-0.14467 


20 








-0.14467 


exact 








-0.14466 
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FIG. 2. Comparison between the Ewald expansion and 
the method of Ref. 18 for Vo4,(r). 



exact values Vo4 a nd Vq g can be used for cubic crystal 
electric field, Eq. ( |2.39a| ,b). For an fee crystal with the 
cubic lattice constant a / 1 we have B4 — \Z4nVo4/a 
and Bq = s/inVoe/a. 

Although we have analyzed both methods at the close 
contact radius R = \/2/4 it is also possible to assess 
the validity of the two methods by decreasing the sphere 
radius r < R. Both methods are designed to describe 
the potential in the interstitial region. They can not 
approximate the potential correctly in the whole space. 
Therefore, both of them are expected to fail at a certain 
small radius, but their departure from the exact solution 
is an indication of their accuracy. The results for Vo4,(r) 
and V e(r) {K max R=90.73, G=40, n=84) at r < R are 
shown in Fig. 2 and 3. 

Weinert' procedure gives more smooth values for r 
close to R, but the Ewald expansion is more stable in 
the whole range. The calculation of the monopole term 
yout ^ Qr j g Q £ S p ec j a i interest. By decreasing the 

sphere radius r we change the density of the compensat- 
ing negative charge distribution. As a result, both Voir) 



FIG. 3. Comparison between the Ewald expansion and 
the method of Ref. 18 for Vo&(r). 



and Vo°"*(r), Eq. ( |3.2| ), are functions of r. From the last 
value however it is possible to reconstruct V^°"* (R) at the 
close contact radius R. On can easily obtain the result 
from electrostatic considerations: 



v out (R) r = 



V out {r) + 2n{R 2 -r 2 )p 
1 - |7r(i? 3 - r 3 )p ' 



where 



1 4tt - 

'7 T r 



(3.7a) 



(3.7b) 



In principle, V^ ut {R) should be independent of r, but in 
practice it exhibit s suc h dependence as indicated on the 
l eft h and side of ( 3.7a ). This is because the estimation 
( |3 . 7a ) depends on the accuracy of the calculated value 
Vq uZ (r) and as such, it is a function of r. We find it 
convenient to study Vg OM *(i?)| r which is expected to be 
constant rather than the initial value Vq iU (r) . In order to 
compare the two techniques we plot the function A (r) = 
{V out \ r -V out \ ex )/V ° ut \ ex mFig. 4. (V out \ ex = -5.5612 
is the best value for Vq u1 calculated from data of Table I, 
last row.) As one clearly sees from Fig. 4, the Ewald 
expansion performs better than the approach of Weinert. 



B. Periodic Coulomb potentials of the I 
cubic charge density 



4 and I 



Here we illustrate how the proposed method works for 
high order multipoles with / = 4 and 1 = 6. We consider 
two cases. In the first case the charge distribution on the 
touching sphere surfaces is given by 



6(r - R) 
P4{r) = 5 K 4\ r h 



(3.8a) 
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FIG. 4. Comparison between the Ewald expansion 
(G = 40, KR < 90.73) and the method of Ref. 18 for V (r) 
(n = 84), see text for details. 



In the second case it is 



6(r - R) 

Pe( r ) = 12 K ^ T > 



(3.8b) 



Here S(r') is the one-dimensionaLdelta-function and K4, 
Kq are cubic harmonics (SAFs)O Now it is not neces- 
sary to introduce a compensating charge in the intersti- 
tial region since t he to tal c harge of the unit cell given 
by integration of (3.8a) or ( |3.8b ) over the polar angles 
yields zero. According to the treatment of Sec. II we 
start with the potential of a single sphere, say S(n = 0). 
The potentials then are 



U 6 s (r) 
for r < R, and 



4tt 1 / r \ 4 

Tr(r) k ^ 

13 R \RJ v ; ' 



Ul(r) 



9 R \r 
4tt 1 /R\ 7 



13 R V r 



(3.9a) 
(3.9b) 

(3.10a) 
(3.10b) 



for r > R. Notice that the potentials Ug and Ug have 
discontinuous derivatives at R. (Here we do not make any 
distinction between U(r) and U(r), Sec. II, since there is 
no renormalization of the charge density inside a spher e.) 
Using the method of Sec. II we find from Eqs. ( 2.21a , b) 
the Fourier transforms of the potentials of the periodic 
arrangements of such charged surfaces on the fee lattice, 



Vi(K) = 
V 6 (K) = 



(4tt) 2 R 
9v K 
(4ir) 2 R , . 



13 v K 



(j 3 (KR)+j 5 (KR))K 4 (K), (3.11a) 
(j 5 (KR)+j 7 (KR))K e (K). (3.11b) 



Notice, that we arrive at the same results if first we 
Fourier transform the densities, 



47T - 

Pi (K) = —j 4 (KR) K A (K), (3.12a) 
v 

p 6 (K) = -—j 6 (KR) K 6 (K), (3.12b) 
and then calculate the potentials by means of 

47T 

MK) = J^P4(K), (3.13a) 

MK) = ^P6(K). (3.13b) 



(In oder to establish the equality with ( 3.11a|,b) we used 
the property 10.1.21 of Bessel functions, Ref. |3l].) As we 
have many times mentioned before (Sec. II) the Fourier 
expansion is valid only in the interstitial region. 
Inside a sphere S(n) the first potential is given by 

Vi(f) - V40 + v M (-^) V 4 (f) + V A6 (-^) V 6 (f) + .... 

(3.14) 

For the second potential, accordingly, we have 

V 6 (r) = V 60 + V 64 (-^) V 4 (f) + V 66 (J)'if,(f) + .... 

(3.15) 

The constants V40, V60, and V44, V46, Vg4, Vg6 are found 



from the boundary conditions, Eqs. (2.24a) and (2.24b), 
where V(K) is given by V 4 (K), Eq. (3.11a), for the first 



potential and by V§{K), Eq. ( 3.111 ), for the secon d po - 
tenti al. (C ompare these results with the potentials (3.9a) 
and (3.9b) for a single sphere.) 

In order to distinguish the potential due to the multi- 
pole moment of the sphere S[n) from the potential of the 
rest of the crystal, we subtract from V4 and V 6 the con- 
tribution i nduce d by the multipole moment of the sphere 
(see Eqs. ( |2.31a -c)) and introduce the parameters 



V'OUt 
44 



v 66 



V44 



66 



We have calculated V 4 °4 Ut , 



4?r 1 

~9R'- 
4tt 1 

13^' 

\rout 

v m 1 



(3.16a) 
(3.16b) 



V46 and Vgo and then 
compared them with those obtained by the method of 
Ref. [l8], see Tables V and VI. Alternatively, we can 
compute and V S"* b y perfo rming the two-center 

multipole expansions ( 2.32 ), ( 2.33 ). To do it, we firs t cal- 
culate integrals W44(n, nf) and V6e(ft,n'), Eq. ( 2.33| ), be- 
tween the functions K4(n) and ^(fi'), or between K§{n) 
and Kq(ji') at sites n and n'. By summing over n' on the 
fee lattice we estimate V44" and V 6 g u *, i.e. 



V'OUt 



v tt (n,n'), 



(3.17) 



where t = 4 or 6. Clearly, the accuracy depends on the 
number of sites n' included in the summations. The re- 
sults of these calculations are given in Table VII. (The 
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representative coordinates of five nearest shells as well as 
the number of atoms of the shells are quoted in Tab. III.) 



TABLE V. Comparison of the direct Fourier expansion 
and the method of Ref. 18. 714 and ne are the parameters 
of Ref. 18 for Va and V&. "Exact" corresponds to Weinert' 
expansion with K max R < 124.08 and 714 = 96, = 94. 





direct expansion 




Ref. 


§ 






Trout 
V 44 


Trout 

v m 


7I4 


Trout 
V 44 


n 6 




35.19 


1.8876 


1.1750 


24 


2.2139 


22 


1.5161 


57.37 


2.0113 


1.3074 


46 


2.2139 


44 


1.5334 


79.60 


2.0688 


1.3696 


69 


2.2140 


67 


1.5220 


101.8 


2.1007 


1.4036 


91 


2.2140 


89 


1.5211 


exact 


2.21402 


1.52112 




2.21402 




1.52112 



TABLE VI. Comparison of the direct Fourier expansion 
and the method of Ref. 18. 714 and n§ are the parameters 
of Ref. 18 for V4 and Vs. "Exact" corresponds to Weinert' 
expansion with K max R < 124.08 and 714 = 96, ne = 94. 





direct expansion 




Ref. 








Vm> ■ 10 2 


Veo 


714 


V ie ■ 10 2 


716 


Veo 


35.19 


-2.5778 


-0.1426 


24 


-1.6667 


22 


-0.1442 


57.37 


-1.9364 


-0.1431 


46 


-1.7042 


44 


-0.1464 


79.60 


-1.7908 


-0.1436 


69 


-1.6820 


67 


-0.1448 


101.8 


-1.7529 


-0.1439 


91 


-1.6812 


89 


-0.14467 


exact 


-1.68119 


-0.14466 




-1.68119 




-0.14466 



From the Tables V-VII it follows that the method of 
Ref. ^| gives a better conv ergenc e despite the fact that 
the expressions ( 3.12a ,b), ( 3.13a ,b) for the Fourier co- 
efficients are exact. Interestingly, the exact expressions 
formally correspond to m = n§ = -1 in Eq. 28 ofRef.pJ, 
while the better convergence is achieved for large positive 
114 and ne implied by the technique of Ref. Q We ascribe 
this property to the fact that the potentials V4 and V@ 
have discontinues derivatives at the sphere boundaries. 
Therefore, the exact expansions (fU2aJ,b), ( pU3a| ,b) try 
to reproduce these cusps while the pseudo-charge den- 
sity of Ref. 18 has n — 1 continuous derivatives thereby 
avoiding t he pr oblem . Thu s, our test calculations for the 
densities ( 3.8a ) and ( 3.8b ) indicate that Weinert' proce- 
dure gives a better convergence. 

We can also generalize the present consideration for 
arbitrary mul tipolc sphere moments Qi(R) and Qe(R) 
given by Eq. (2.13a) on an fee lattice when its cubic lat- 
tice constant a / 1. For the former case we obtain that 
the potential inside a sphere due to the other spheres' 
multipoles Q4 of the fee crystal is 



r \ 



V: ut (r) = B i0 + B i4 [-\ # 4 (f) + fl 46 - K 6 (r) 



RJ 



RJ 



(3.18) 



where -B40 = k^y^, B44 = k^V^ 1 , B 4e = &4V46, and the 
coefficient of proportionality for touching spheres (R = 
v^a/4) is given by 



64 

k A = -zQiiR). 

n 



(3.19) 



Analogously, the expression for the potential Vg ut due to 
Qe(R) is obtained as 



TABLE VII. Calculation of V44 and Vee by summation 
over neighboring shells. The discrepancy between the con- 
verged value for Vm and the exact one from Table VI is due 
to numerical errors of the two center integration. 



shell 


7J44 


Trout 
V44 


«66 


Trout 
V 66 


1 


0.17897 


2.14767 


0.12632 


1.51593 


2 


0.01406 


2.23204 


0.00064 


1.51976 


3 


-0.00091 


2.21024 


0.00005 


1.52088 


4 


0.00035 


2.21444 


0.00002 


1.52106 


5 


-0.00002 


2.21393 


0.0 


1.52102 


10 




2.21403 




1.52102 


20 




2.21401 




1.52102 


exact 




2.21402 




1.52112 



R 



V 6 ou \r)=B 60 + B M {-) K4f)+B 66 {-) K 6 (r 



\R 



(3.20) 



where B 60 = k e V<w, B 64 = k e V 6 4, B e6 = /c 6 V 6 ° 6 u *, and 



^12 

h = —Qe(R). 



(3.21) 



Notice that here we can immediately determine the po- 
tentials using exact values V40 = V04, V44, V45, Vgo = ^06, 
Vq4 and Vqq quoted in Tables V and VI. These parame- 
ters as well as those calculated in Sec. III. A are in fact 
important structural constants of the fee lattice. 



IV. HARTREE-FOCK-ROOTHAAN METHOD: 
CALCULATIONS OF THE DIRECT MATRIX 
ELEMENTS 



The Haptf|ee-Fock operator for the electronic system is 
defined asl3'El 
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T = E (Hd + Jd~ ICd) , 



(4.1) 



where Ti.d is the Hamiltonian operator for a dth electron 
moving in the field of nuclei alone, 



7~Ld — —- 



V 



E 



t \Rd ~ X(n) 



Jd and ICd are the direct Coulomb and exchange opera- 
tors, respectively. In order to obtain the best trial func- 
tion of the dth electron, ipd{x) = {x\d), one considers 
basis functions Xa( x ) = (x\a), an d expands ipd in terms 
of Xa, 



4>d{x) 



Ida Xa{x), 



(4.3) 



where 7^ are the coefficients of the expansion, and the 
index x stands for the coordinates R and the spin pro- 
jection s z = ±1/2 of the electron, x = (R,s z ). For the 
matrix elements of Jd and ICd one has 

{a\Jd\b) [ X*a(xi)Xb(xi)^ip*(x 2 )'ipc{x2)dx 1 dx2, 



(4.4a) 



(a\fC d \b) = y~] / x*a{xi)i'c{xi)-^—Xb{x2W c {x2)dx 1 dx2- 

(4.4b) 



In Eqs. (4.4a,b) the summation is understood over all 



occupied electron states except d. In the following we 
consider—the restricted HF method. Usually in the HF 
methodtll the wave function i^ d ( x ) — <?K-R)C( s z) is given 
by the coordinate (orbital) part cj>(R) and the spinor part 
C(s z ). However, it is well known that the spin-orbit 
coupling mixes the two components £ + = £(1/2) and 
£_ = £(—1/2). (The spin-orbit coupling is especially im- 
portant for core shells.) Therefore, here we consider more 
general spin-orbitals, 

M x ) = MR)(+ + <f>-(R)(- (4-5) 

Using the time reversal symmetry one can construct 
from ijid the time reversed state tp d = Ttpd (see Appendix 
C). According to the Kramers theorem ip d and i\)* d have 
the same energy and are orthogonal. Therefore, such 
treatment is analogous to the conventional restricted HF 
methodtll with the double occupancy of e^. However, in 
the following we will treat these two states separately, 
as different components of a doubled valued irreducible 
representation. Therefore, although we work within the 
restricted HF met hod , there are no factor 2 in front of 
Jid a nd Jd in Eq. (|4.l|). Notice, that as a consequence of 
(4.5), the integrals ( E4aj ,b) include summations over two 
spin components. 
We regroup T as 



t = r 



COUI rp&XC 



where 



(4.2) and 



T Coul = Ypi d + J d ), 



(4.6) 
(4.7) 

(4.8) 



The matrix elements of J :Coul can be written in the fol- 
lowing form: 

(a\T Coul \b) = J X * a (x) (-^V 2 + V d (R^j X b(x)dx, (4.9) 

where the Hartree (electrostatic) potential V d acting on 
the dth electron reads 



v\r) = y' [ r ^6 x>) dx> y 

~ J \R - R n ^ 



z 



t \R-X(n) 



It is convenient to rewrite the potential as 
V d (R) = V(R) - V' d (R), 



(4.10) 



(4.11) 



where V(R) is the potential of all electrons and nuclei 
considered in section II, 



V(R) = 



p(R') 

\R-R' 



-dR' 



(4.12) 



while V' is the potential created by the dth electron 
alone, 



J \R-R'\ 
The charge density p(R) in ( 4.12| ), 
P (R) = p el {R) + p n (R), 
comprises the point charges of nuclei, 

n 

and the total electronic charge density, 

p el (R) = j2r c (x)Mx)- 



(4.13) 



(4.14a) 



(4.14b) 



(4.14c) 



(Again, in (4.14c) the summation is implied over two spin 
components.) In solids we are dealing with the two types 
of electron states, which are extended states of valence 
electrons (with the total charge density p va i(R)) and lo- 
calized states of core electrons (p core (f(n))), and 
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Pel 



(R) = Pva i(R) + Y,Pcore(R~ X(n)). (4.14d) 



D: 



(4.17) 



The Hartree-Fock-RoothaaH, method implies a self- 
consistent field procedure. In particular, for itinerant 
states (d — (fc,a)) one solves the secular equation 

^2 J r ba {k,a)^ a (k,a) = E(k, a) ^ Oba(k) 7 (fc, a). 

a a 

(4.15) 

Here k is the wave vector, a is the band index, O is 
the overlap matrix and E(k, a) is the corresponding HF 
energy.EZl The HF operator T depends implicitly on the 
solution (via the coefficients j a ). Below in this section 
we calculate the matrix elements of J? Coul j while in Sec. 
V we consider the matrix elements of J- exc . 

A. Multipole expansion of electron density 



In order to calculate the matrix elements of !F 



Coul 



for 

a dth electron from Eqs. (4.9) we have to kn ow the elec- 
trostatic potential V d , Eq. ( ]4 lOp or ( |4.1l| ). We have 
already thoroughly studied this problem in Sec. II and 
know how to proceed starting with the dual representa- 
tion (2.2a,b) of charge density. The only problem is to 



construct the charge density of the itinerant and the lo- 
calized electrons. In this subsection we do it explicitly 
for the core (Al) and the valence (A2) electrons. 

Al. The density of the inner closed shells at site n 
reads 



where 



Pcore{r{n)) = y"\i r p T (f(n)), 



p T (^ = \(f,s z \r)\ 2 = \(x\r)\' 



(4.16a) 



(4.16b) 



Here (x\t) stands for the wave functions of localized elec- 
trons, and n T = 1 if e( r) < E p and zero otherwise (Ep is 
the Fermi energy). In ( }4.16a ) summation over s z — ±1/2 
is implied. Here we will not consider the case of partially 
filled core shells. (Such situation occurs in lanthanides 
with localized 4/ electrons or in transition elements with 
d electrons.) Furthermore, we assume that all core elec- 
trons at site n are confined inside the sphere S(n) so that 
the core wave functions with their first derivatives fall to 
zero at R. 

The inner closed shells are classified according to the 
principal quantum number n, the total angular momen- 
tum (including spin) j and the orbital number I. In the 
case of the spherical symmetry these (2j + 1) electronic 
states belong to a doubled valued irreducible represep, 
tation Dj of the 3-dimensional rotation group SO(3)£B 
In the presence of crystal environment these levels (ex- 
cept s) in general are split into doubled valued irreducible 
representations (r, v) of a crystal double group]23 



where v is used to label representations which occur more 
than once. We shall classify these core electronic states 
according to the atomic indices n, I, j, and T, v 1 k (k 
labels the rows of T) i.e. r = (n, T, v, k), and 

(x\t) =K T (r){r,s z \T) = TZ T (r){f,s z \n,l,j; r,v,k). (4.18) 

The orientational (spin-orbital) functions are linear com- 
binations of real SAFs and spinors, 

(f\n,l,j; r,^,fc> = ^c r (A, S2 )S A (r!)C( S ,). (4.19) 

In the case of the spherical symmetry A = (I, m), r = 
(n,j,nij) aad-c T (A, s z ) are Clebsch-Gordan (or Wigner) 
cocfficicntsEjnj but for other point groups these coeffi- 
cients are not well known. We have derived c T (A, s z ) 
for the cubic double group 0' h in Appendix C. Func- 
tions (x\t) = (f, s z \t) are independent of the radial part 
Ttn,j,i{r), which is assumed to be the same for all states 
belonging to the same n, j and I. The electronic states 
| to, I, j; r, v, k) distinguished by the index k have the same 
energy, e(l,j;T, u) — e(r). We nex t consider the multi- 
pole expansion of p T (f), Eq. (1.16b), and find 



P 



(4.20) 



Since the core density is invariant of the point group, only 
the full symmetrical irreducible representations has to be 
considered, with T\ — A\ g . The coefficients ca(t, t') are 
given by 

+ 1/2 . 

ca(t,t')= J2 {n,s z \T)*S A (0.){n,s z \T')da (4.21) 

s z = -l/2 J 

(We recall that £1 or f stand for the two polar angles 
(9,0).) 

A2. The density of the valence electrons reads 



Pval 



(4.22) 



where (x\ka) is the wave function of a delocalized elec- 
tron with the wave vector k and the band index a: nr 
is the occupation number, 

We expand the electron density of itinerant electrons 
in multipole series inside the spheres. In the following we 
will use the LAPW basis functions, Appendix B. Using 
Eqs. (jB8| ) we find that the local density at a site n is 
given by 



Pval 



(f(n)) 



A 



PA(r(n))S A (f(n)), 



(4.23) 



where 
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(4.24) 



Here and below summation is understood over the re- 
peated indices p and t, and 



Pa - N 2^/ n k,a ( 



(h,p)(h,t) 



(k,a;k,a). (4.25) 



(See Appendix B for definitions.) The coefficients 



A 



(k, a; k, a) are given by 



lh,p)(h,t) 



(k,a; k,a)= 



Ai(ii)A 2 (i 2 ) 



x[ 7 A]P;(fc,a, Sz )[ 7 A] t A2 (fc,a, S ,) C A(A 1 , A 2 ). (4.26) 

We recall that A = (l,T, v, k) (Appendix A) and the 
summations in (4.26) are performed over the subindiccs 
r, v, k of A(/i) and A(Z 2 ) within the manifolds l\ and l 2 . 
Finally, 



c A (Ai, a 2 ) = / s Xl (n)s A (n)s X2 (n)dn. 



(4.27) 



These coefficients can be tabulated before the self- 
consistent-ficld HFR procedure. The density of conduc- 
tion electrons stays invariant under all symmetry oper- 
ations, which means that in expansion (4.27) we con- 



sider only irreducible representations of A\ g symmetry, 
Ta = Ai g . In such case the nonzero coefficients in (4.27) 
can occur only if (1) T\ 1 = T\ 2 or (2) both T\ 1 and T\ 2 
belong to the A\ g irreducible representation. 

In the interstitial region p va i (R) is expanded in Fourier 
series, 



Pval (R) = Pval (K) e 
K 



iKR 



(4.28) 



where from Eqs. (Bl) and (B7) we find 

k,a K' s " 

X7^,(fc, a, s z ) lK+K>( k , s z ). 



(4.29) 



B. Direct Coulomb matrix elements of core states 

If one wants to use the Hartree-Fock-Roothaan method 
then the problem of the radial dependence of core states 
arises. In LDA it is solved quite naturally since it is pos- 
sible to introduce a self-consistent spherically symmet- 
ric potential which includes an average exchange term. 
Nonspherical contributions are small and usually omit- 
ted. Then the radial components are obtained through 
the solution of the Schrodinger (or Dirac) equation in 
the potential. However, there is no such convenience in 



the Hartree-Fock approach since the exchange contribu- 
tion generally can not be reduced to an effective single 
particle potential. In the HF method one routinely uses 
radial dependencies of Gaussian or Slater-type (GT or 
ST). An alternative choice of complete radial basis func- 
tions is given in Appendix D. In all cases the radial part 
is approximated as 



(4.30) 



where u n is a radial basis function (GT, ST or other), 
while r\ stands for the parameters specifying the function. 

The orientational vectors (x\j, l,rrij) are usually spher- 
ical harmonic spinorsB a In a crystal, the spherical sym- 
metry is reduced, Eq. ( 4.17 ), and we replace the spinors 
by their symmetry adapted combinations (x\t) as we dis- 
cussed in Sec. IV Al and Appendix C. Then the basis 
function of a core state r reads 



X T n {x) = (x\t,t)) = u T n {r)(x\T). 



(4.31) 



Since the orientational vector (x\t) is independent of 77, 
the basis functions are distinguished only by the radial 
component u^(r). 

In order to calculate the matrix elements of the di- 
rect Coulomb interaction for each atomic orbital r we in- 
troduce the Hartree potential so that no self-interaction 
occurs, 

V T (r) = V(R) - V T (r) = VJ (r) + V& (r) , (4.32) 

where VJ and refer to the spherical and the nonspher- 
ical components of V T , respectively. First we consider the 
spherically symmetric component, 



VJ(r) = V (R) - VXf). 



(4.33) 



We get 



(T,r,\rf oul \T, V >) = (t, V \T\t, V >) 



+ / u T Jr)uL(r)V ( ;(r)r 2 dr. (4.34) 



Here (r, t]\T\t', 17') are the kinetic energy integrals which 
are well known for GT (ST) orbitals. For the alterna- 
t ive s et of basis radial functions (Appendix D) instead of 
(4.34) the matrix elements are given by Eq. (D3). 



In order to compute the matrix elements of VL we 



partition V(R) in two parts, Eq. fl2.25| ). For Us (the 
potential of the charges inside the sphere S(n)) we dis- 
tinguish further two con tributi ons, from the itinerant and 
the core electrons, Eq. ( 1.14d ). Using equations derived 
in Sec. II, we obtain 

(T, V \FSr l \T, V ') = J2'cA(r,r)(Dr'(out) 

A#0 

+ DY , '(val) + D™'(core)), (4.35) 
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where 



DY> (out) = 



irout 

A 7)7) 



Ri 



-1i 



(4.36) 



Here V£ ut is given by (2.28b) and the multipolar mo- 
ments </j are 



cj," I ul(r)ul,(r)r l+2 dr. 



(4.37) 



For the contribution from the extended states, 
{vat), we obtain 

^»"'M)=^Q(«|«,)P^' P)(i2 '* ) ! (4-38) 



hh 



where the two- fold radial integrals C/(...) are defined by 
@ and p^'P^ 2 -^ is defined by Eqs. ( pL2"i| , ( gjg ). The 



contribution from the other core electrons reads 

Q(uX'I^T'ftr') ca(t, r') <,. (4.39) 



D^' 77 (core 



E ( 



The last sum runs over all occupied states except r (no 
self-interaction), i.e. n T r = and ri£ = 1 if the core state 
t' 7^ r is occupied and zero otherwise. 

In cubic crystals only / and d shells are split by CEF 
effects, Appendix C. The CEF interaction of lo calize d 
electrons with the delocalized ones (second part of_( 4.35])) 
has been discussed in a number of papers, Refs. H-Elll 
and calculated in Refs. |36|,|37]. The important result here 
is that we have obtained all CEF interactions. In partic- 
ular, besides those considered in Refs. |3^,^, [WCrijiclude 
the CEF effects from the rest of the crysta£3~H3 (first 
contribution, Eq. ( [4.35 )) and from the o ther core shells 
with nonspherical density (third part of ( 4.35| ) ) . 



C. Direct Coulomb matrix elements of extended 
states 



Here we consider matrix elements of F Coul ; Eq. (4.7), 
for a conduction electron \d) = \k, a). Usually in the 
LAPW method the basis functions k(R) = {R\k,K) 
(see Appendix B) are defined in an effective potential 



LDA 







Vn 



-Vq xc which includes the direct Coulomb,ip. 



teraction Vo and the LDA exchange potential l/ ea:c .Eirt£l 
The potential V^ DA is spherically symmetric inside the 
"muffin-tin" (MT) spheres and is a constant in the in- 
terstitial region. Notice that V$ DA is used only for the 



construction of the basis functions Xk k(R), Eq. (Bl). 
Having defined x% jt{R)> m principle one can calculate 
the matrix elements of a general potential V tot (R) = 
V Coul (R) + V^LRlJthe procedure is known as FP- 
LAPW methodESO'Ea). Below we follow the same ap- 
proach, but instead of Vq DA for construction of the ba- 
sis functions we consider only the electrostatic potential 
Vq(R) without any exchange, 



V Q d (R) = V a (R)-V d (R). 



(4.40) 



Here V' (R) is the spherically symmetric Coulomb poten- 
tial due to the electron d. As we will see later in section V 
we can omit V' {R) for any conduction electron, because 
the corresponding matrix elements, Eq. ( p. 48 ), decreases 
as 



(k,K'\V d \k,K)~^ 



(4.41) 



and vanish in the limit N — > oo. Therefore, constructing 
the LAPW basis functions x £ it > we can use ^ ne potential 
V (R), i.e. Vq(R) = V (R). Next step is to calculate the 
direct matrix elements for the conduction electron d in 
the full Coulomb potential V(R). Following the method 
described in Sec. II we write V(R) = V {R) + V M {R) 
for the potential inside the spheres (Vq and Vm are the 
spherically symmetric part and the contribution due to 
the other multipoles, correspondingly). In the intersti- 
tial region we use the Fourier expansion of V{R). Then 
the Hartree-Fock operator j: Coul is separated into three 
parts, 



r 



Coul 



rr-Coul 

— • / " 



^Coul 
f M 



1~Coul 
' •'I > 



(4.42) 



where T§ oul comprises the kinetic energy and the spher- 
ically symmetric potential Vq(R), F^° ul accounts for the 
other multipole terms while ^ oul stands for the elec- 
trostatic interaction in the interstitial region. Starting 
with Tq ou1 we arrive at the standard expressions for the 
matrix elements of the LAPW methodjlj 

(p, K', s' z \^ oul + Kj\k, K, s z ) = Hg,#(js) (4.43) 

where Hg, g(k) is given by Eq. (15) of R ef. |l5| . (The 
overlap matrix is given by Eq. (13a) of Ref. jig .) 

The other contributions of the general Coulomb po- 
tential V Coul (R) follow from equations of Sec. II. After 
some algebra, we obtain: 

(p, K', s' z \T%r l + F? ou % K, s z ) = 6j: 6 SzS , (4.44) 



>A#0 hh 



{h,p){h,t)-(h,p){h,t) + Bi 



Here Bj arises due to the interstitial contribution, 
Bi = Y; V (P) °(K - K' + P), 



(4.45) 



where V(P) stands for the Fourier coefficients, Eq. 
(|2.22b|) . 0{P') is given by © and 



Y, E A x * i (k,K')A{ 2 {k,K)c A (X 1 ,\ 2 ). 

Ai(ii) A 2 (/a) 

(4.46) 
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To condense notations we shall use L for (l,p). In (4.44) 



Vd,(R) 



b l ± l 2 = b^ iL2 (out) + B^ lL2 (val) + B^ lL2 (core), (4.47) 

where the index out stands for the contribution from the 
other spheres, 



Pcb(R') 



dR' 



-trout 



(4.48) 



\R-R>\ 

which corresponds to the "exchange" density 

1/2 

Pcb(R)= ^2 i/j*(R s z)xb(R,s z ). 

s»=-l/2 



(5.2) 



(5.3) 



Here V^" is given by (2.28b) and the multipolar mo- 
ments lL2 are 



(■0 C has two spin components due to the spin-orbit inter- 
action, Eq. (4.5).); 2) We calculate the matrix element of 
exchange as 



1l 



(h, P )(h,t) 



(4.49) 



(a\T^\b) =J2j V cb {R)p* ca {R)dR. 



(5.4) 



For the contribution from the delocalized electrons we 
obtain 



(4.50) 



where p L ^ L ' 2 and Cj(...) are given by Eqs. ( 4.25 ) and (El), 
respectively. 

For the contribution from the closed shell core elec- 
trons we have 



We want to stress that the "exchange" potential V c b and 
the "exchange" density p c b are technical quantities here, 
which are employed only for t he c alculation of the ma- 
trix element of exchange, Eq. (5.4). They should not be 
confused with the effective exchange potential which is 
derived and widely used in the local density approxima- 
tion. 



B 



(h,p)(h,t) 

A 



(core) =J2 C li u h u U n rRr)cA(r,T)n T , (4.51) 



wher e ag ain the two-fold radial integral Ci is defined by 
Eq. (||l[). 

The non-s pheri cal components of the total Coulomb 
potential in ( 4.44 ) represent crystal field like effects for 
conduction electrons. 



V. CALCULATIONS OF THE EXCHANGE 
MATRIX ELEMENTS 

In this section we calculate the exchange matrix ele- 
ments of the Hartree-Fock, method. The most general 
expression for exchange isO 



A. Exchange for localized electrons 

Here we calculate the exchange for a localized state 
r sited at ft, i.e. d = r. Index c can refer either to a 
conduction state or to another core state r'. First we 
consider c as a conduction state, i.e. |c) = \k, a). The ra- 
dial wave function of the l ocalized electron r is expanded 
in terms of uL, Eq. ( 4.30 ), so that \a) and \b) stand for 
|t, 77). Since the state r is co nfin ed inside S(n), for the 
calculation of the exchange ( |5.4| ) we need to know the 
density p c b and the potential V cfc , Eqs. ( f>.3\ ) and (pT^), 
only inside the same sphere. For p cb we get 



1 



P(k,*)(T, v )(r(fi)) = t=^Pa >(«))SA(r(r?)), 



(5.5) 



{a\T exc \b) 



X*a{x 1 )i3 c {xi)^-Xb{x 2 )lp* c {X2)dx 1 dX2- 

(5.1) 



where 



(k,a) (r,f)) / 



Here both indices a and b refer to basis wave functions of 
an electronic state d (whic h ca n be either a conduction 
or a localized state), Eq. ( |4.3|) ; ip c stand for the esti- 
mated wave functions of conduction and localized elec- 
trons obtained from a previous iteration of the HFR self- 
consistent procedure. The summation is understood over 
all occupied stat es c. As before (Sec. Ill), x = (R, s z ) and 
the integration (5.1) includes summation over two spin 
components s z . We will calculate the matrix element 
fl5.1| ) in two steps: 1) we construct an auxiliary Coulomb 
potential 



Pa' "'i r ) = 5Z£ CA ( r ' Sz;A ') h A ]\*(^ a ^ s z) 

x«?(rK(r). (5.6) 
The coefficients ca(t, s z ; A') are given by 



CA (r, Sz - a) = / (si, Sz \T)s A {n)s x {n)dn. 



(5.7) 



From (5.5) the exchange between |r) and \k,a) is found 
as 

(T,r,\^ xc (k,a)\T, V ') = ^D^{k,a), (5.8) 
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where 

D™ c (k,a)=J2Y. E ca(t, S2 ;A 1 ) Ca (t,4;A 2 ) 

A Ai A2 s z ,s' z 

x[ 7 A]l*(k,a,s z ) [ 7 ^' 2 (fc,a,^)G(«|<<0- (5-9) 



Here the integrals C;(...) are given by Eq. (|El[). The 
exchange with all extended states reads 

(r, V \^ XC (c)\r,v') = ^5>~ c & (5.10) 



Next we consider the exchange (^J) between r (as be- 
fore, r = d) and a core electron r' localized at the same 
site n. r' is now our reference state c, i.e. c = r'. Pro- 
ceeding analogously, we find that the multipole expansion 
of the "exchange" density p c b is 

Pr' { r,r,){r{n)) = £ P A ' ^ (r(r?)) 5 A (r(n)), (5.11) 

A 

where 

^' (T "' ) (r)=c A (r ; rO^(r)<(r). (5.12) 
The exchange integral is 

(r, r?| ^ ( CO re)|r, rf) = ]T <'^ c (t'), (5.13a) 

t' 

where 

£ exc (r') =^|cA(r,r')| 2 Q(^^|^ r m;,). (5.13b) 



B. Exchange between an extended state with the 
localized states 



In this subsection we consider d as an extended elec- 
tron state, i.e. d = (k, a), while c refers to a core state r 
localized inside a sphere S'(n), c = r(n). We expand 
the wave function in terms of (R, s z \k, K) = 

(R\k,K)((s z ), Appendix B. Therefore, a,b — (k,K,s z ). 
The exchange density p c b, Eq. (|5.3|), is located inside the 
sphere S(n), 

ViV A A' 

xc A (T, Sz ;X')A p x ,(k,K)) S A (f), (5.14) 



where the coefficients ca(t, s z ; A') are given by Eq. (p.7[). 
Again, using t he m ultipole expansion ( |5 . 14 ) we calculate 
the exchange (5.4) and obtain 



where 

B exc (r) 



A Ai A 2 



c A (t, s*; Ai) ca(t, s' z ; A2) 



x a»; (fc, k)4 2 (fc, /?') c, (7e T < |k t < ) . (5.16) 

Notice that the latte r res ult is independent of n and 
the matrix elements (5.15) can be nonzero even for off- 
diagonal spin functions (s z = —1/2, s z = 1/2 or vice 
versa) due to the spin-orbit interaction. In order to ob- 
tain the exchange with all localized electrons, we sum 
( |5.15 ) over all sites n and all occupied core electron states 
|t), and find that 

(k,K, Sz \T exc (core)\k,K' lS ' z )=J2B exc (T). (5.17) 



C. Exchange between extended states 

The calculation of exchange between two delocalized 
electrons is much more i nvol ved and quite laborious. As 
the state c in Eqs. ( 5. 1 )-( 5.4 ) we consider now a conduc- 
tion state \p, (3) with the wave function ip c {R) = (x\p, /?). 
As before, d = (fc, a) and the wave function of the dth 
electron, (R\k, a), is expanded in terms of LAPW ba- 
sis functions (R, s z \k, K) — (R\k, K)((s z ) (Appendix B) 
labeled by indices a,b — (k,K,s z ). In this subsection 
we first calculate the matrix elements of dth electron ex- 
changed with the extended state c. Then by summing 
over all occupied extended states |c) = \p, (3) we will be 
able to compute the matrix element of exchange between 
d and the other conduction electrons. 

Inside a sphere S(n) the "exchange" density p c b, Eq. 
T3I), is given by 



Pc(k,K,s z ) 



(f) 



1 



A 

(5.18) 



where q=p—k and 

p f,K,s s){r) = J2J2 ca(Ai,A 2 ) l-yA]l*(p,f3,s z ) 

I1I2 A1A2 

xA{ 2 (k,K)ul(r)ul(r). (5.19) 
Outside the spheres p c b is expanded in plane waves: 
1 



Pc(k,K,s z ) 



[R] 



Nv 



s z e 



iK'-R 



(5.20) 



K' 



Proceeding as in Sec. II we continue the plane wave 
representation ( |5.20|) i nside the spheres and then subtract 
it out from Eq. (|5.1q). This procedure renormalizes the 



(k,K,s z \^ c (r(n))\k,K', S ' z ) = -B^(r), (5.15) are given by 



multipole radial functions p c ^- k,K,Sz ^ f ( 5.19| ), which now 
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p'^ k ' g '"\r) =p c A {k ' K ' s *\r) 



v 



K' 



Here S\(K' +q) is understood as the symmetry adapted 
function S&_(Clg, + g) of the polar angles defined by 

the vector K 1 + q. 

Notice that p c rg g s ) is not a periodic function of R. 
A translation by t = X(n) transforms it as 



{ Pc(k,K,sJ R ) 



(5.22a) 



This transformation is not identical for q ^ 0. From Eq. 
fl5.2| ) we find that the same transformational law holds 
for the corresponding potential, i.e. 



c(k,K,..)&) = e imH) K(k,K,sjR)- (5.22b) 



As a consequence of Bloch's theorem we obtain 



(k,K,s z ) 
*c(k,K,s z ) 



(R) 
(R) 



e m P f t a 

" c \k,K,- 



— e iqR 



*c(k,K,s z ) 



(R), 



where p c b and V c b are periodic functions, i.e. 



(5.23a) 
(5.23b) 

(5.24a) 
(5.24b) 



Therefore, in the interstitial region p c b and V c b are ex- 
panded in Fourier series, while the initial functions p c b 
and V c b are expanded in terms of exp(i(q + K')R), where 
K' is a reciprocal lattice vector. The factor exp(i<fi?) or, 
more precisely, exp(igA), will be also present for solu- 
tions inside the spheres. Indeed, if one knows a solu- 
tion inside a s phere, say S{il = 0) at the origin, then by 
means of Eq. (5.22b) it is easy to generate the solution 
inside any other sphere. Notice, however, that the factor 



eyLj>(iqX (ri)) cancels in the final expression (5.4). As a 
result, it is easily to generalize the method of Sec. II for 
the present consideration. 

First we consider the potential in the interstitial region, 



y e (i2) = le^5>e(i?')< 



AK'R 



N 



(5.25) 



K' 



As in Sec. II we distinguish three contributions there, 
V e (K') = V*(K') + V e s '°(K') + V e s > M {K'), (5.26) 



where Vf'°(K') and Vf' M (K') are the Fourier compo- 
nents of the potentials of the monopole and the other 
higher multipoles of the sphere, respectively. Vj(K') rep- 
resents a component from the the plane wave expansion 
(^20|), i.e. 



Vj(K') = - 



4tt 



v\K' + q\ 



1%-,, if(P,P,Sz)- 



2 'K'+K 



(5.27) 



The Fourier components Vf<°(K') and V e s ' M (K') are 
found as 

V e S ' 9 (K') = - f e- l ^ +fl "> n Vs^R)d 3 R, (5.28) 

v J cell 

where g = 0, M, and 



Here lAg a and are the potentials due to the monopole 
(/ = 0) a nd th e other multipoles (/ > 1) of a sing le sphere 
(see Eq. (|2~10|)). For the potential^, Eq. ( |2~To| ), we use 
the equivalen t cha rge distribution given by Eq. ( |5.21 ). 
The integral ( 5.28| ) is taken over a unit cell of the Bravais 
lattice. One can show that Vf' 9 (K') is also expressed 
asS 



V e s ' 9 (K') = ~ f e-^+*' )n U Sg {R)d 3 R, (5.30) 

^ J crystal 

where the integration spans the whole crystal. The rep- 
resentation ( |5.30 ) is then used to compute the Ewald 
expansion coefficients V e s -°(K') and to find V e S ' M (K') by 
direct integration. The procedure is the same as in Sec. 
II. 

In order to proceed with the Ewald expansion, we in- 
troduce an effective point "exchange charge" inside the 
spheres, 



4?r / p 





,c(k,K,s z ) 



(r)r 2 



(5.31) 



Using the orthogonality of the radial functions u\ and ■ 



and the relation co(A, A') = <5xa'/v47t, we obtain 



-vC (k.K ,s z ) 

to 



I X 



4ttR 2 ^ j^K* + g]R) , 

Is l£»; ■ * 7^, + ^(P. (5-32) 



'• \ K '+<i\ 



where A'; 1 = 1 (normalization of the radial functions of 
the LAPW method) and 



Nf = / [uf(r)} 2 r 2 dr. 



(5.33) 



Although in general QQ^ k ' K,s *' (R) ^ in Appendix F we 
prove that for the present particular case (one atom in the 
unit cell) this "charge" has to be considered only for q = 
0, i.e. when p = q but a ^ (3. However, the final results 

will not be affected if one assumes that Ql {k ' K ' s *\R) ^ 
for q 7^ 0. If there are two or more atoms in the unit 
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cell, the "exchange charges" of nonequivalent atoms must 
be introduced and taken into account even for the case 
p^q, Appendix F. We shall now proceed further as for 
the general case. The Ewald expansion coefficients are 
given by 



where U ei s(r(n)) is the potential of the charges located 
inside the sphere S(n) and V°^{r) accounts for the po- 
tential of the rest of the crysta l. Th e single sphere po- 
tential U ev s(r) is given by Eq. ( 2.26 ) where qo(r), q' (r) 
are replaced by 



D -\K'+q\ 2 /AG 2 n c(k,K,Sz) 
Wo 



(5.34) 



« \K' + q\ 2 

For the higher multipoles we introduce two functions 



,c(k,K,s z ), s n 

q Q '(r) = V47r 



^c(k,K,s z ), , I ,c(k,K,s z ), r ,1+2 , , 

h ( r ) = P A {r)r dr 

Jo 

y c(k,K,s z ) , \ / / c(k,K,s z ) I i \ ,1—1 j i 



(5.35a) 
dr. (5.35b) 



Proceeding then as in section II we find: 

V A#0 

X ( 41 ^ ^ { ) \K' + q\R^ 



' pf^^r'dr', 

J r 

and q A (r), <?A'( r ) b Y 
Jo 

f R _ - 

,c(k,K,s z ) f s _ / c(fc,R>*)^ ,1 

9 a v) — Pa 



V)' 



(5.41a) 
(5.41b) 

(5.41c) 
(5.41d) 



(5.36) 



The potential of the "exchange charges" out side t he 
sphere 5(n) is found by subtracting from Eq. ( 5.39a ,b) 
the potential of the "charges" inside. Thus, we introduce 
the constants 



where 



AA(\K' + q\) 



<:.. 



c(k,K,s z ) 



(R) 



V e °X = Ve,A 



R 

4tt q c A &K' s '\R) 
21 + 1 &+ 1 



(5.42a) 



(5.42b) 



xji(\K' + q\r)r 2 dr. (5.37) 
Eqs. ( pTTj ), fl5.34|) and (|5.36| ) fully determine the plane 



wave components V e (K'), Eq. ( |5.26[ ). Thus, we have ob- 
tained the potential V e (Rj, Eq. ( |5.25| ), in the interstitial 
region. 

The potential inside a sphere S(n) is found as 
V e (r) = ie^W [V e ,o(r) + ^V e , A (r) S A (f)], (5.38) 

A#0 

where V^oM and V EtA {r) are smooth functions of r. As 
before, we introduce the constants V e ,o = V e ,o(R) and 
V e> A = V e! A(^)- They are found from the boundary- 
value problem for the surface of the sphere: 

V e , =J2io{\K' + q\R)Ve(K'), (5.39a) 

K' 

Ve, A = 4m l Y,3l(\K' + q\R) S A (K> + q) V e (K'). 

(5.39b) 



Then the potential inside the sphere due to the "exchange 
charges" outside is 



A#0 



(5.43) 



A"' 



Here V C: o and V &: a depend on q and have to be recalcu- 
lated for each extended state p — k + q. We recall that 
Ve(K') is given by the sum ( |5.26| ). Following Sec. II. C 
we rewrite the potential V e (r), Eq. (5.38|), as 



Therefore, Eqs. ( |5~4"ol ), ( ^43| ) and (|2^ ) where q a (r), 
q' (r), qA{r), 9A'( r ) are given by ( |5.41a -d) fully determine 
the effective "exchange" potential inside any sphere. 

In order to obtain the matrix element of exchange we 
integrate the "exchange" potential with the "exchange" 
density p* t - (r), Eqs. (E>.18])-(|5.2C|), and recall that c 

stands for (p,(3). From the potentials inside the spheres 
and in interstices we distinguish three contributions to 
the exchange, 

(k,P,s' z \T exc {p^)\k,K^ z ) = l(BZmP) 

+Br c {p 1 P) + Bf c (j> 1 p)). (5.44) 

Here B^, Bg xc and Bf cc are the contributions from in- 
tegrations with V°g t (r), U e ,s(r) inside a sphere, and with 

V e (R) in the interstices, respectively. First we carry out 
the integrations inside a sphere S(n), and then perform 
summation over the N spheres. As a result we get 



Ve{r{n)) = -e™ [U e , s (r(n)) + V°$(f{n))] , (5.40) 



A#0 



q C A (kM (R) 
R l 



(5.45) 
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The contribution from the interstitial region gives 
BT C (P,P) = E V e (K') 7p+it (p,0,s' z )O(K> - P), (5.46) 

where 0(K") is the overlap integral ( |3.5| ). Finally, for 

Bg xc we get 

C(^) = EEEE E CA(A 1 ,A 2 )c A (A' 1) Ay 

x [ 7 A)l* (p, (3, s' z ) [ 7 A}(, (p, f3, s z ) A{ 2 (k, K) a{,* (k, P) 

xQ(ulul\u(4 2 ). (5.47) 

Notice that all contributions to the Fock exchange, Eqs. 
(p4|)-(p6|) and ( f)A7\ ) are proportional to 1/N as one 
could exp ect f rom a charge distributed over N cells. 
From Eq. ( 5.44 ) we can calculate the Coulomb self-energy 
associated with a conduction state \k,a). Assuming 
p = k, a = J 3 and taking into account the expansion 
(B6) and (B7), we arrive at 



E si (k, a) = V h + E E 7 K fc > a ' ^i<( k : a ' 



K,P t 



x {k,P,s' z \F™ c {k,a)\k,K,s z ) ~ -, (5. 



where (k, P, s' z \T exc {k, a)\k, K, s z ) is given by ( 5.44 ). 
Here we have introduced Vh which is an electrostatic en- 
ergy, associated with the electron charge e distributed 
homogeneously inside the crystal. (Such homogeneous 
charge distribution is absent for exchange if two con- 
duction states are different, Appendix F. It is compen- 
sated by the positive contribution of nuclei for the direct 
Coulomb interaction.) Vh depends on the shape of a crys- 
tal, but decreases as As a result we observe that 
E s i — > in the limit N — > 00 for any extended state. 

By summing the exchange over all occupied conduction 
states \p, (3) ^ \k, a) we obtain a finite value of exchange 
for each conduction electron |fc,a), 

(k, P, s' z \T exc (val)\k, K, s z ) = 1 J2( B out(P> P) 

+Br c (p,(3) + Br(p,(3))- (5.49) 

In the absen ce of t he sp in-orbit coupling, the matrix ele- 
ments ( [5.44 ) and ( [5.49| ) are diagonal in spin components 
s z = ±1/2. An important practical complication here 
arises due to the fact that the structural constants V e o, 
V e \ (VgQ*, V°\^) depend on p (q) and the calculation of 
them has to be repeated for each vector p. 



VI. CONCLUSIONS 

We have presented a new Hartree-Fock-LAPW method 
for electron band structure calculations. The method 



combines the restricted Hartree-Fock-Roothaan ap- 
proach with the crystalline basis functions in the form 
of linear augmented plane waves.. The strategy of the 
full potential LAPW treatmenttSEfl is adopted for cal- 
culations of the matrix elements of the direct Coulomb 
interactions and exchange. This is pivotal for collecting 
all exchange terms together including the long-range and 
multipole contributions. 

In the framework of the FP-LAPW treatment an orig- 
inal technique for the solution of periodic Poisson's equa- 
tion is formulated, Sec. II. The technique takes into ac- 
count the partitioning of space into two regions, inside 
the spheres and in the interstices, Fig. 1. In the intersti- 
tial region we expand electron densities and the poten- 
tial in Fourier series and express "exchange" densities in 
terms of plane waves. Inside the spheres we expand den- 
sities and potentials in multipole series. Finally, we use 
these expansions to calculate the matrix elements of the 
direct Coulomb interaction (Sec. IV) and the exchange 
(Sec. V). The crystal field effects are considered for core 
electron shells and for conduction electrons. These effects 
are associated with the nonspherical density components 
of /, d and, for noncubic symmetries, of p electrons £3cil 
There, the crystal site symmetry is taken into account 
and the basis functions are adapted for the spin-orbit 
interaction. 

The technique for solving Poisson's equation has been 
applied to the face centered cubic lattice, Sec. III. We 
have calculated structural constants which are used to 
restore cubic Coulomb potentials inside a sphere from its 
monopole (/ = 0) and multipole (I = 4, 6) moments. We 
have compared our technique with the pseudo-charge- 
density method of Weinert, Ref. [l8] and the two-center 
expansion of the Coulomb interaction, Ref. 

At present we are working on programming the formu- 
las derived in this article. However, it is already clear 
that the task consists of two independent parts. First 
of all, one should calculate the multipole matri x ele- 
ments ca(Ai,A2) of electron transitions, Eq. ( 4.27 ), and 
the other coefficients related to them (such as ca(t, t'), 
Eq. ( |4.2l| ), and c A (r, s z \ A), Eq. (|J)). The integrations 
there involve only angular (and spin) parts of electronic 
wave functions and thus the coefficients can be tabu- 
lated and stored before the HFR self-consistent proce- 
dure. Also to this part one should add calculations of 
the relevant structural constants, such as V40, V^"*, V46, 
^60, and V 6 g ut computed in Sec. Ill for the face centered 
cubic structure. (The constants are needed to restore 
the full potential for a given set of multipole moments.) 
These calculations depend on the type of crystal symme- 
try but are separated from the problems of HFR method. 
The second task is to program the matrix elements and 
all relevant procedures of the restricted HFR method. 
For those purposes one can start with an existing code 
of LAPW method and develop it on the basis of the con- 
siderations presented in this article. 
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APPENDIX A: 

Throughout this paper instead of complex (surface) 
spherical pharmonics YJ m we use real symmetry adapted 
functionscEl (SAFs) 5a, which transform according to ir- 
reducible representations A(Z) of a site symmetry group 
G. The composite index A(i) stands for (l;T(l), k) 
where T(l) labels the irreducible representations within 
the I manifold, fj,(l) numbers the representations that oc- 
cur more than once and k denotes the rows of a given 
representation. In the following we omit index I in T and 
A. SAFs Sa are linear combinations of Y™ with the same 
I, i.e, 



i 

m=0 



/ 



(Al) 



m— 1 



where the coefficients c„ c and c m k 



depend on the group 
"' "•' - Y°. The coef- 



G under consideration, and Y l 
ficients c s for different groups are quoted in Tables 
2.4-2.6 of Ref. ||. There are (21 + 1) independent SAFs 
5a belonging to the I manifold. The real spherical har- 
monics are 



Y 



Y 



1 

i 

iV2 



(Y t m + Yr m ) 



y; 



(A2a) 
(A2b) 



where F ; m are taken with the phase definition of Ref. |2^. 
(It is different from the definition used by Condon and 
Shortly.Ej) For a given I in a {21+1) dimensional space we 



consider row vectors (Y, ' c ,Y l 



.,Y^) 



0,c 1A l,c 

Y (we exclude Y^^ = 0), and (S\) = S, and the matrix 
= c. Then the SAFs and the spherical harmonics 



for a given I are connected through an orthogonal trans- 
formation, 

S = Y-c. (A3) 
One can easily find the inverse transformation since 



Y = S-. 



(A4) 



where T stands for the transpose, since c _1 = c T . It is 
more convenient to use Sa than Y t m due to their known 
symmetry properties. 



In the density expansion, Eq. ( |2.2b| ), only the SAFs 
of A\ g symmetry survive because density stays invari- 
ant under all symmetry operations of G. However, for 
calculations of exchange (Sec. V) there is no such simpli- 
fication and the full basis set (including SAFs belonging 
to the other irreducible representations) should be taken 
into account. 

We use SAFs to describe both electronic densities and 
wave functions. For electronic states we adopt a nota- 
tion with small letters, i.e A = (l;T, v, k). For localized 
electrons for conciseness we incorporate also the princi- 
pal quantum number n and write r = (n, Z; T, u, k). In 
the latter case T refers to a double valued irreducible 
representations of G.Ej 



APPENDIX B: 

Here we introduce some definitions and notations pf, 
the linear augmented plane wave method (L AP W) .E£nL£l 
The coordinate basis functions are plane waves in the 
interstitial region, 



x% AR) = {R\k,K) 



J(k+K)-R 



'Nv 



(Bl) 



and a linear combination of local atomic functions inside 
the spheres, 



(R e S\k,K) 



n A 



(B2) 



wherel 



§ 



^P(f) = (A x {k,K) Ul {r) + Bxik.K)^)) S x (f). 

(B3) 

(If one uses spherical harmonics instead of SAFs S\, 
then A = (l,m).) In order to condense notations we 
introduce two components (p=l,2) of the radial function 

uf(r), 

uUr)=Mr), uUr)^u l (r)= dU f E E) - (B4) 

and the corresponding to them two components 
AR(k,K), which are 

A\ {k,K) = A x (k,K), A{ {k,K) = B x {k,K). (B5) 

The coefficients A p x are obtained by requiring that the 
basis functions and their derivatives are continuous on 
the sphere boundary.lij 

In the absence of a static magnetic field and the spin- 
orbit coupling, the conduction electronic states with spin 
projections s z = ±1/2 are degenerate. (The spin-orbit 
coupling can be included later in the second variational 
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treatment as described in Ref. |l(|) The wave function 
of a conduction electron juitii the wave vector k and the 
band index a then readsacia 

1/2 

(R\k,a)= J2 (R,s z \k,a)C(s z ), (B6) 

s z = -l/2 

where £(s 2 = ±1/2) are the two spinors and 

(R,s z \k,a) = ^ 1& {k,a,s z ){R\k,K). (B7) 
it 

The coefficients ^(fc, a, s z ) are found by the HF varia- 
tional procedure. E3 Inside a sphere S(n) the wave func- 
tion is given by 

(r,s z \k,a) = j= ^2[yA} p x (k,a,s z )uf(r) S x (f), 



(B8) 



N 



where r = r(n) and summation over p = 1,2 is implied. 
Here we have introduced the notation 

[lA] p x (k, a, s z ) = 1k$> «> s z )A p x (k, K). (B9) 
it 



the first, Eq. (C2a), through the time reversal symmetry. 
One can easily check it by noting that the SAFs are real 
and by applying the following rules for the time reversal 
symmetry 



For two components of F 8 we have 

(r\d 5 /2, (r 8 , 1)) = ^=(i25 2i (T2 g ,3) + 3S 2 .(E g .2))(+ 
V 15 



(C3) 



+ ^TF( _i5, 2.(T2 g .l) + S'2,(T2 g ,2))C-, (C4a) 

V 15 



(f|d 5 /2,(r 8l 2)) 



V5 



'5'2,(_Eg,l)C+ 

(^2,(T2 ff ,l) +5 , 2,(T2g,2))C-- (C4b) 



The other two components are obtained from ( C4a| ,b) by 
employing the rules (C3). 

We consider analogously the splitting of the f 5 j 2 (I = 
3) states in the cubic symmetry and obtain 

1 i— 2 

(^1/5/2, (r 7 , 1)) = ~j={-ViiS 3 ,.A2u + ^-5*3,(T2 M ,3))C+ 

2 

"I ^=(5 , 3,(T2u,l) + ^>3,(T2ii,2))C-) (C5) 

V 21 
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and 



Here we derive analytical orientational wave vectors 
for the cubic site symmetry Oh by employing the eigen- 
vectors tabulated in Ref. |3q . 

The core states of S1/2, P1/2, P3/2 an d ^3/2 electrons re- 
main degenerate while d 5 / 2 , fs/2 and / 7 / 2 are split in the 
cubic environment. If -D5/2 and D7/2 are doubled valued 
representations of 50(3), then the symmetry lowering 



5/2 

V/2 



r 7 + r 8 , 
r 7 + r 8 + r 6 . 



(Cla) 
(Clb) 



For two components of the doublet r 7 of d 5 / 2 we have 
found 



1 

(f|d 5/2 ,(r 7 ,l)) = -j=S 2 ,(T2g,3)(+ 



+ -^(iS 2 ,(T2g,l) ~ S 2 ,(T2g,2))C-> (C2a) 



(f\d 5/2 , (r 7 ,2)) = -j={iS 2 ,(T2g,l) +52,(T2 fll 2))CH 



(C2b) 



where 5 l 2j (T2g,fe) refer to the three components (k = 1 — 
3) of T2 9 symmetry (7 = 2) giv en in Table 2.6 of Ref. 
p9| . The second component, Eq. (C2b), is connected with 



/ 5 3 

(r\h/2, (r 8 , 1)) = -y — S , 3,(T2«,3)C+ + [^-^(^.(Tlu.l) 

1 / 5 

— *S3,(TlTt,2)) + 2 V 21 ^ 3 .(T2«,1) + ^3,(T2u,2))]C-) (C6a) 

(rl/5/2, (r 8 ,2)) = y / ^s , 3i(T i^ 3) c+ - ^=[\/3(5 3i(T i Ui i) 

(C6b) 



) - \/5(S 3 ,(T2m,1) - iS 3 ^T2u,2))]C- 



Rcal SAFs 5 3 {k = 1 - 3) 

are given in Table 2.6 of Ref. ^9|. The other components 
are obtained through the time reversal symmetry, Eq. 



Finally, for / 7 / 2 core states we get 

1/7/2) ( r 6, 1)) = -7=(£3,(Tlu,l) - «S , 3,(Tlu,2))C+ 



1 



^^3,(ri«,3)C-) 



(C7a) 



(»* 1/7/2) ( r 7, 1)) = ~^{'2.iS 3 ,A2u + S3,(T2«,3))C+ 

+ "-yy(5 , 3,(T2ii,l) + «S' 3 ^T2tt,2))C-- (C7b) 

These are first components of Tg and T 7 , respectively. 
For the quartet T 8 of / 7 / 2 we have found 
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1 / 5 

(r|jV/2, (r 8) 1)) = 21 ^3,(T1«,1) - iS'3,(Tl«,2)) 
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3 /"g - The two-fold integral C;(/i, /a 1 5102) of four radial 

"^(^3,(T2t»,i) + *<S , 3,(T2m,2))]C+ + y 2j^3,(Tiu,3)C-) (C8a) f unc tions /i(r), / 2 (r), .9i( r ) an d 52(0 is denned as 



(f|/ 7/2 ,(r 8 ,2)) 



7 -5*3,(T2u,3)C+ + gtV ^(^3,(Tlu,l) G (/i , / 2 |0i5 2 ) 



-^3,(Tlu,2)) — \/ y (5 , 3,(T2«,1) + *<S3,(T2u,2))]C- 



(C8b) 



<ir r 2 / cZr' r' 
Jo 

x/i(r)/ 2 (r)« ; (^^).9i(^').92(r') (El) 



Again, the other components of Tq, IV and Tg, ar e fo und 
by ap plying the time reversal transformation (C3) to 
dC7jb) and (ggaLb). 



where the one center multipole function is given by 
vi(r,r') = -^i^tt, Eq. (2.£), and r< (r>) is the smaller 
(larger) of r and r' . 
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APPENDIX F: 



Here we describe a simple method to generate radial 
basis functions. We can use the spherically symmetric 
component of the total electrostatic potential Vo(r) (in- 
side a sphere S(n)) to find the radial solution u n j and 
the corresponding energy E n j, 



hiu n j — EnjUnj. 



(Dl) 



Here n is the principal quantum number, and hi is the 
radial operator of the Schrodinger equation (Eq. (lb) of 
Rcf. [j^). We assume that the solutions are confined in- 
side the sphere S(n) and on the sphere boundary r = R 
we have 



ui(R) = 0, 



dui(r) 



Or 



0. 



(D2) 



The boundary co nditio ns a rc complementary to (Dl). 
Starting with (Dl) and (D2) one obtains u n ,i and E n< i. 

For a given / we consider the radial functions u n> i which 
differ from each other by the principal quantum number 
n = I + 1, 1 + 2, ... . These functions correspond to the 
same angular dependence (specified by I, mi) and form a 
complete orthonormalized set. Therefore, they can serv e 
as a basis for any radial function 1Z T satisfying Eq. (|I5^) . 
A nonrelativistic function of a localized electron is char- 
acterized by the combined index r = (I, mi) and we write 
uZ = U nt i, i.e. rj = n. For the relativistic case one has 
to distinguish two solutions with the total momentum 
j = 1-1/2 and j = 1 + 1/2. Then r = mj) and the 
basis functions are u^, where again rj = n. 

For the matrix eleme nts of the spherically symmetric 
potential Vq , Eq. (4.33), we find 



77|^ Go "'|r,7/) = E n .j6 w - — C «X^ T ft T ). (D3) 



Here E n> j is the eigenvalue corresponding to vL (J 
(l,j)) and the integral Cq(...) is given by Eq. (Ell). 



We consider the exchange between two extended states 
\d) = \k,a) and \c) — \p,0). We exp and both states in 
the LAPW basis functions, Eq. (B7). Proceeding as in 
section V.C, we introduce an effective "exchange" den- 
sity p c d(R), Eq. (5.3), and the corresponding "exchange 
charge" : 

Qo d (K) = E EEW(P- ft s s ) hA]l(k, a, s z ) Nf 



+0(p,/3;k, a), 



where 
and 



v ^ \K' + q\ 



(Fl) 

pg,(p,P;k,a) 
(F2) 



p Rl {p,l3]k,a) = ^y^2,'y*g., + j t {p,l3,s z )^g{k,a,s z ). (F3) 
R s ; 

On the other hand, the orthogonality relation for the two 
extended states is 

(p,0\k,a) = S aP S fi : 

l X(l) 

+0(p, (3; k, a) + 5^ £ ]T 7 * g (p, (3, s z ) lR {k, a, s z ). (F4) 
it s - 



Comparing it with (Fl) we observe that if p ^ q then 
)Q d (R) — 0. However, it will not be erroneous to use 



this charges as they appear in HFR procedure, Eq. (5.32) 



Since Poisson's equation is linear, their contributions will 
cancel in the fin al results. Generalizing the orthonormal- 
ity relation (F4) for the case of few atoms one can show 
that it leads to 
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£0? 



0, Py^k. 



(F5) 



We conclude that the "exchange" charges Q\ d are not 
necessarily zero. (Here i labels different atoms in the unit 
cell.) The orthogonality relation ( |F5| ) ensures that there 
is no Coulomb divergence associated with the uniform 
component of the "exchange density" . 
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